!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

module cfc11_mod

!BOP
! !MODULE: cfc11_mod
!
! !DESCRIPTION:
!
! !REVISION HISTORY:
!  SVN:$Id: cfc11_mod.F90 28439 2011-05-18 21:40:58Z njn01 $

! !USES:

   use POP_KindsMod
   use POP_ErrorMod
   use POP_IOUnitsMod

   use blocks, only: nx_block, ny_block, block
   use domain_size, only: max_blocks_clinic, km, nx_global, ny_global
   use domain, only: nblocks_clinic, distrb_clinic
   use exit_mod, only: sigAbort, exit_POP
   use communicate, only: my_task, master_task
   use prognostic, only: tracer_field
   use kinds_mod
   use constants, only: c0, c1, char_blank, delim_fmt
   use io, only: data_set
   use io_types, only: stdout, datafile, io_field_desc, io_dim,       &
       nml_in, nml_filename, construct_file, construct_io_dim,        &
       construct_io_field, rec_type_dbl, destroy_file,                &
       destroy_io_field, get_unit, release_unit
   use io_tools, only: document
   use tavg, only: define_tavg_field, accumulate_tavg_field
   use timers, only: get_timer
   use passive_tracer_tools, only: forcing_monthly_every_ts,          &
       init_forcing_monthly_every_ts, ind_name_pair, tracer_read,     &
       rest_read_tracer_block, file_read_tracer_block, read_field

   implicit none

  !-----------------------------------------------------------------------------
  !   public/private declarations
  !-----------------------------------------------------------------------------

   private

! !PUBLIC MEMBER FUNCTIONS:

   public :: &
       cfc11_tracer_cnt, &
       cfc11_init, &
       cfc11_set_sflux,  &
       cfc11_tavg_forcing, &
       SCHMIDT_CFC

!EOP
!BOC
!maltrud
     character(char_len) :: mmessage

  !-----------------------------------------------------------------------------
  !   module variables required by forcing_passive_tracer
  !-----------------------------------------------------------------------------

  integer(int_kind), parameter :: &
       cfc11_tracer_cnt = 1

!-----------------------------------------------------------------------
!     index bounds of passive tracer module variables in TRACER
!-----------------------------------------------------------------------

   integer (int_kind) ::  &
     num_cfc11_years

  !-----------------------------------------------------------------------------
  !   relative tracer indices
  !-----------------------------------------------------------------------------

  integer(int_kind), parameter :: &
       cfc11_ind = 1     ! cfc11 index

  !-----------------------------------------------------------------------------
  !   derived type & parameter for tracer index lookup
  !-----------------------------------------------------------------------------

  type(ind_name_pair), dimension(cfc11_tracer_cnt) :: &
       ind_name_table = (/ ind_name_pair(cfc11_ind, 'CFC11') /)

  !-----------------------------------------------------------------------------
  !   derived type for tracer initialization
  !-----------------------------------------------------------------------------

    type(tracer_read) :: &
         gas_flux_fice,       & ! ice fraction for gas fluxes
         gas_flux_ws,         & ! wind speed for gas fluxes
         gas_flux_ap            ! atmospheric pressure for gas fluxes

  type(forcing_monthly_every_ts), save :: &
       fice_file,              & ! ice fraction, if read from file
       xkw_file,               & ! a * wind-speed ** 2, if read from file
       ap_file                   ! atmoshperic pressure, if read from file

  !-----------------------------------------------------------------------------
  !   define tavg id for 2d fields related to surface fluxes
  !-----------------------------------------------------------------------------

   integer (int_kind) :: &
      tavg_FICE,         &! tavg id for ice fraction
      tavg_XKW,          &! tavg id for xkw
      tavg_ATM_PRESS,    &! tavg id for atmospheric pressure
      tavg_pCFC11,       &! tavg id for cfc11 partial pressure
      tavg_PV,           &! tavg id for piston velocity
      tavg_SCHMIDT_CFC11,&! tavg id for cfc11 schmidt number
      tavg_CFC11SAT       ! tavg id for cfc11 saturation

  !-----------------------------------------------------------------------------
  !   define array for holding flux-related quantities that need to be time-averaged
  !   this is necessary since the forcing routines are called before tavg flags
  !-----------------------------------------------------------------------------

   real (r8), dimension(:,:,:,:), allocatable ::   &
      CFC11_SFLUX_TAVG

   integer (int_kind) :: cfc11_sflux_timer

   character(char_len) :: &
         cfc11_formulation,           & ! how to calculate flux (ocmip or bulk)
         cfc11_interp_weight_file,    & ! filename for hemispheric weights
         cfc11_interp_weight_file_fmt,& ! file format for weights file
         cfc11_time_series_file         ! filename for ascii time series of atm cfc11

  logical(log_kind), dimension(:,:,:), allocatable, save :: &
       LAND_MASK

   real (r8), dimension(:,:,:), allocatable :: &
      CFC11_INTERP_WEIGHT,  &! latitudinally dependent weight for atm CFC11
      pCFC11                 ! surface patial pressure of CFC11

   real (r8), dimension(:,:), allocatable :: &
      cfc11_time_series ! time series of atm CFC11

!EOC
!***********************************************************************

contains

!***********************************************************************
!BOP
! !IROUTINE: cfc11_init
! !INTERFACE:

 subroutine cfc11_init(init_ts_file_fmt, read_restart_filename, &
                       tracer_d_module, TRACER_MODULE, errorCode)

! !DESCRIPTION:
!  Initialize cfc11 tracer module. This involves setting metadata, reading
!  the module namelist and setting initial conditions.

! !REVISION HISTORY:
!  same as module

! !USES:

   use broadcast, only: broadcast_scalar, broadcast_array
   use constants, only: c0, field_loc_center, blank_fmt, &
       field_type_scalar
   use prognostic, only: nx_global, ny_global, curtime, oldtime
   use grid, only: KMT, zt, zw, n_topo_smooth, fill_points
   use forcing_tools, only: find_forcing_times
   use time_management, only: freq_opt_nyear, freq_opt_nmonth

! !INPUT PARAMETERS:

   character (*), intent(in) ::  &
      init_ts_file_fmt,    & ! format (bin or nc) for input file
      read_restart_filename  ! file name for restart file

! !INPUT/OUTPUT PARAMETERS:

   type (tracer_field), dimension(cfc11_tracer_cnt), intent(inout) :: &
      tracer_d_module   ! descriptors for each tracer

   real(r8), dimension(nx_block,ny_block,km,cfc11_tracer_cnt,3,max_blocks_clinic), &
      intent(inout) :: TRACER_MODULE

! !OUTPUT PARAMETERS:

   integer (POP_i4), intent(out) :: &
      errorCode         ! returned error code

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

    character(*), parameter :: subname = 'cfc11_mod:cfc11_init'

    character(char_len) :: &
         init_cfc11_option,           & ! option for initialization of bgc
         init_cfc11_init_file,        & ! filename for option 'file'
         init_cfc11_init_file_fmt       ! file format for option 'file'

    integer(int_kind) :: &
         n,                   & ! index for looping over tracers
         k,                   & ! index for looping over depth levels
         l,                   & ! index for looping over time levels
         ind,                 & ! tracer index for tracer name from namelist
         iblock,              & ! index for looping over blocks
         nml_error              ! namelist i/o error flag

    type(tracer_read), dimension(cfc11_tracer_cnt) :: &
         tracer_init_int,      &! namelist variable for initializing tracers
         tracer_init_ext        ! namelist variable for initializing tracers

    namelist /cfc11_nml/ &
         init_cfc11_option, init_cfc11_init_file, tracer_init_ext, &
         init_cfc11_init_file_fmt, cfc11_interp_weight_file,       &
         cfc11_interp_weight_file_fmt, cfc11_time_series_file,     &
         gas_flux_fice, gas_flux_ws, gas_flux_ap, &
         cfc11_formulation

   character (char_len) ::  &
      cfc11_restart_filename  ! modified file name for restart file

!-----------------------------------------------------------------------
!  initialize forcing_monthly_every_ts variables
!-----------------------------------------------------------------------

   errorCode = POP_Success

   call init_forcing_monthly_every_ts(fice_file)
   call init_forcing_monthly_every_ts(xkw_file)
   call init_forcing_monthly_every_ts(ap_file)

!-----------------------------------------------------------------------
!  initialize tracer_d values
!-----------------------------------------------------------------------

   tracer_d_module(cfc11_ind)%short_name = 'CFC11'
   tracer_d_module(cfc11_ind)%long_name  = 'CFC11'
   tracer_d_module(cfc11_ind)%units      = 'fmol/cm^3'
   tracer_d_module(cfc11_ind)%tend_units = 'fmol/cm^3/s'
   tracer_d_module(cfc11_ind)%flux_units = 'fmol/cm^2/s'

    !---------------------------------------------------------------------------
    !   default namelist settings
    !---------------------------------------------------------------------------

    init_cfc11_option = 'unknown'
    init_cfc11_init_file = 'unknown'
    init_cfc11_init_file_fmt = 'bin'
    cfc11_interp_weight_file = 'unknown'
    cfc11_interp_weight_file_fmt = 'bin'
    cfc11_time_series_file = 'unknown'
    cfc11_formulation = 'unknown'

    do n = 1,cfc11_tracer_cnt
       tracer_init_ext(n)%mod_varname  = 'unknown'
       tracer_init_ext(n)%filename     = 'unknown'
       tracer_init_ext(n)%file_varname = 'unknown'
       tracer_init_ext(n)%scale_factor = c1
       tracer_init_ext(n)%default_val  = c0
       tracer_init_ext(n)%file_fmt     = 'bin'
    end do

    gas_flux_fice%filename     = 'unknown'
    gas_flux_fice%file_varname = 'FICE'
    gas_flux_fice%scale_factor = c1
    gas_flux_fice%default_val  = c0
    gas_flux_fice%file_fmt     = 'bin'

    gas_flux_ws%filename     = 'unknown'
    gas_flux_ws%file_varname = 'XKW'
    gas_flux_ws%scale_factor = c1
    gas_flux_ws%default_val  = c0
    gas_flux_ws%file_fmt     = 'bin'

    gas_flux_ap%filename     = 'unknown'
    gas_flux_ap%file_varname = 'P'
    gas_flux_ap%scale_factor = c1
    gas_flux_ap%default_val  = c0
    gas_flux_ap%file_fmt     = 'bin'

    if (my_task == master_task) then
       open (nml_in, file=nml_filename, status='old',iostat=nml_error)
       if (nml_error /= 0) then
          nml_error = -1
       else
          nml_error =  1
       endif
       do while (nml_error > 0)
          read(nml_in, nml=cfc11_nml,iostat=nml_error)
       end do
       if (nml_error == 0) close(nml_in)
    endif

    call broadcast_scalar(nml_error, master_task)
    if (nml_error /= 0) then
       call document(subname, 'cfc11_nml not found')
       mmessage = 'ERROR : stopping in '/&
                                         &/ subname
       call exit_POP(sigAbort,mmessage)
    end if

    !---------------------------------------------------------------------------
    !   broadcast all namelist variables
    !---------------------------------------------------------------------------

    call broadcast_scalar(init_cfc11_option, master_task)
    call broadcast_scalar(init_cfc11_init_file, master_task)
    call broadcast_scalar(init_cfc11_init_file_fmt, master_task)
    call broadcast_scalar(cfc11_interp_weight_file, master_task)
    call broadcast_scalar(cfc11_interp_weight_file_fmt, master_task)
    call broadcast_scalar(cfc11_time_series_file, master_task)
    call broadcast_scalar(cfc11_formulation, master_task)

    call broadcast_scalar(gas_flux_fice%filename, master_task)
    call broadcast_scalar(gas_flux_fice%file_varname, master_task)
    call broadcast_scalar(gas_flux_fice%scale_factor, master_task)
    call broadcast_scalar(gas_flux_fice%default_val, master_task)
    call broadcast_scalar(gas_flux_fice%file_fmt, master_task)

    fice_file%input = gas_flux_fice

    call broadcast_scalar(gas_flux_ws%filename, master_task)
    call broadcast_scalar(gas_flux_ws%file_varname, master_task)
    call broadcast_scalar(gas_flux_ws%scale_factor, master_task)
    call broadcast_scalar(gas_flux_ws%default_val, master_task)
    call broadcast_scalar(gas_flux_ws%file_fmt, master_task)

    xkw_file%input = gas_flux_ws

    call broadcast_scalar(gas_flux_ap%filename, master_task)
    call broadcast_scalar(gas_flux_ap%file_varname, master_task)
    call broadcast_scalar(gas_flux_ap%scale_factor, master_task)
    call broadcast_scalar(gas_flux_ap%default_val, master_task)
    call broadcast_scalar(gas_flux_ap%file_fmt, master_task)

    ap_file%input = gas_flux_ap

    do n = 1,cfc11_tracer_cnt
       call broadcast_scalar(tracer_init_ext(n)%mod_varname, master_task)
       call broadcast_scalar(tracer_init_ext(n)%filename, master_task)
       call broadcast_scalar(tracer_init_ext(n)%file_varname, master_task)
       call broadcast_scalar(tracer_init_ext(n)%scale_factor, master_task)
       call broadcast_scalar(tracer_init_ext(n)%default_val, master_task)
       call broadcast_scalar(tracer_init_ext(n)%file_fmt, master_task)
    end do

    !---------------------------------------------------------------------------
    !   initialize tracers
    !---------------------------------------------------------------------------

    select case (init_cfc11_option)

    case ('ccsm_startup', 'zero', 'ccsm_startup_spunup')
      TRACER_MODULE = c0
      if (my_task == master_task) then
          write(stdout,delim_fmt)
          write(stdout,*) ' Initial 3-d CFC11 set to all zeros'
          write(stdout,delim_fmt)
          call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
      endif

    case ('restart', 'ccsm_continue', 'ccsm_branch', 'ccsm_hybrid' )

       cfc11_restart_filename = char_blank

       if (init_cfc11_init_file == 'same_as_TS') then
          if (read_restart_filename == 'undefined') then
             call exit_POP(sigAbort, &
                'no restart file to read cfc11 from')
          end if
          cfc11_restart_filename = read_restart_filename
          init_cfc11_init_file_fmt = init_ts_file_fmt

       else  ! do not read from TS restart file

          cfc11_restart_filename = trim(init_cfc11_init_file)

       end if

       call rest_read_tracer_block(init_cfc11_init_file_fmt, &
                                   cfc11_restart_filename,   &
                                   tracer_d_module,          &
                                   TRACER_MODULE)

    case ('file')
       call document(subname, 'cfc11 being read from separate file')

       call file_read_tracer_block(init_cfc11_init_file_fmt, &
                                   init_cfc11_init_file,     &
                                   tracer_d_module,          &
                                   ind_name_table,           &
                                   tracer_init_ext,          &
                                   TRACER_MODULE)

       if (n_topo_smooth > 0) then
          do k=1,km
             call fill_points(k,TRACER_MODULE(:,:,k,1,curtime,:), &
                              errorCode)

             if (errorCode /= POP_Success) then
                call POP_ErrorSet(errorCode, &
                   'cfc11_init: error in fill_points')
                return
             endif
          enddo
       endif

    case default
       call document(subname, 'init_cfc11_option = ', init_cfc11_option)
!      call exit_POP('ERROR: stopping in ' // subname)

    end select

    !---------------------------------------------------------------------------
    !   apply land mask to tracers
    !---------------------------------------------------------------------------

    do iblock=1,nblocks_clinic
    do n = 1,cfc11_tracer_cnt
       do k = 1,km
          where (k > KMT(:,:,iblock))
             TRACER_MODULE(:,:,k,n,curtime,iblock) = c0
             TRACER_MODULE(:,:,k,n,oldtime,iblock) = c0
          end where
       end do
    end do
    enddo

    !---------------------------------------------------------------------------
    !   allocate and initialize CFC11 partial pressure array
    !---------------------------------------------------------------------------

    allocate( pCFC11(nx_block,ny_block,max_blocks_clinic) )
    pCFC11 = c0

    !---------------------------------------------------------------------------
    !   allocate and initialize LAND_MASK (true for ocean points)
    !---------------------------------------------------------------------------

    allocate( LAND_MASK(nx_block,ny_block,max_blocks_clinic) )
    LAND_MASK = merge(.true., .false., KMT > 0)

    call get_timer(cfc11_sflux_timer, 'CFC11_SFLUX',1, &
                                          distrb_clinic%nprocs)

!-------------------------------------------------------------------------------
!  call other initialization subroutines
!-------------------------------------------------------------------------------

    call init_cfc11_tavg
    call cfc11_init_sflux

!-----------------------------------------------------------------------
!EOC

 end subroutine cfc11_init

!***********************************************************************
!BOP
! !IROUTINE: init_cfc11_tavg
! !INTERFACE:

 subroutine init_cfc11_tavg

! !DESCRIPTION:
!  Define tavg fields not automatically handled by the base model.

! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      var_cnt             ! how many tavg variables are defined

!-----------------------------------------------------------------------

   var_cnt = 0

   call define_tavg_field(tavg_FICE,'CFC11_FICE',2,               &
                          long_name='CFC11 Ice Fraction',         &
                          units='fraction', grid_loc='2111')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_XKW,'CFC11_XKW',2,                 &
                          long_name='CFC11 XKW',                  &
                          units='m/s', grid_loc='2111')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_pCFC11,'pCFC11',2,                 &
                          long_name='CFC11 partial pressure',     &
                          units='atmospheres', grid_loc='2111')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_PV,'CFC11_PV',2,                   &
                          long_name='CFC11 piston velocity',      &
                          units='m/s', grid_loc='2111')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_ATM_PRESS,'CFC11_ATM_PRESS',2,     &
                          long_name='CFC11 Atmospheric Pressure', &
                          units='atmospheres', grid_loc='2111')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_SCHMIDT_CFC11,'CFC11_SCHMIDT',2,   &
                          long_name='CFC11 Schmidt Number',       &
                          units='none', grid_loc='2111')
   var_cnt = var_cnt+1

   call define_tavg_field(tavg_CFC11SAT,'CFC11_SAT',2,            &
                          long_name='CFC11 Saturation',           &
                          units='mmol/m^3', grid_loc='2111')
   var_cnt = var_cnt+1

!-----------------------------------------------------------------------

   allocate(CFC11_SFLUX_TAVG(nx_block,ny_block,var_cnt,max_blocks_clinic))

!-----------------------------------------------------------------------
!EOC

 end subroutine init_cfc11_tavg

!***********************************************************************
!BOP
! !IROUTINE: cfc11_init_sflux
! !INTERFACE:

 subroutine cfc11_init_sflux

   use broadcast, only: broadcast_scalar, broadcast_array
   use constants, only: field_loc_center, blank_fmt, field_type_scalar
   use grid, only: KMT, zt, zw
   use forcing_tools, only: find_forcing_times

   type (datafile) ::    &
      in_file             ! data file type for init ts file

   type (io_field_desc) :: &
      io_tracer           ! io field descriptors for input Tracer

   type (io_dim) :: &
      i_dim, j_dim, k_dim, month_dim ! dimension descriptors

   real (r8), dimension(:,:,:,:), allocatable :: &
      TEMP_DATA           ! temp array for reading Tracer data

    integer(int_kind) :: &
         io_error,            & ! io error status
         nu,                  & ! io unit number
         n,                   & ! index for looping over tracers
         k,                   & ! index for looping over levels
         nc,                  & ! index for looping over columns in time series file
         ny,                  & ! index for looping over years in time series file
         iblock                 ! index for looping over blocks

    !---------------------------------------------------------------------------
    !   read gas flux forcing (if required)
    !   otherwise, use values passed in
    !---------------------------------------------------------------------------

   select case (cfc11_formulation)

     case ('ocmip')

       allocate(TEMP_DATA(nx_block,ny_block,12,max_blocks_clinic))

    !---------------------------------------------------------------------------
    !   first, read ice file
    !---------------------------------------------------------------------------

       allocate(fice_file%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
       in_file = construct_file(fice_file%input%file_fmt,             &
                         full_name=trim(fice_file%input%filename),    &
                         record_length = rec_type_dbl,                &
                         recl_words=nx_global*ny_global)
       call data_set(in_file,'open_read')

       i_dim = construct_io_dim('i',nx_global)
       j_dim = construct_io_dim('j',ny_global)
       month_dim = construct_io_dim('month',12)

       io_tracer = &
           construct_io_field(trim(fice_file%input%file_varname), &
           dim1=i_dim, dim2=j_dim, dim3=month_dim,                &
           field_loc = field_loc_center,                          &
           field_type = field_type_scalar,                        &
           d3d_array=TEMP_DATA)

       call data_set(in_file,'define',io_tracer)
       call data_set(in_file,'read'  ,io_tracer)
       call destroy_io_field(io_tracer)

!untested       !$OMP PARALLEL DO PRIVATE(iblock, n)
       do iblock=1,nblocks_clinic
       do n=1,12
          fice_file%DATA(:,:,iblock,1,n) = &
             TEMP_DATA(:,:,n,iblock)*fice_file%input%scale_factor
          where (.not. LAND_MASK(:,:,iblock))
             fice_file%DATA(:,:,iblock,1,n) = c0
          end where
       end do
       end do
!untested       !$OMP END PARALLEL DO

       call data_set(in_file,'close')
       call destroy_file(in_file)

       if (my_task == master_task) then
          write(stdout,blank_fmt)
          write(stdout,'(a12,a)') ' file read: ', &
             trim(fice_file%input%filename)
       endif

       call find_forcing_times(fice_file%data_time, &
            fice_file%data_inc, fice_file%interp_type, &
            fice_file%data_next, fice_file%data_time_min_loc, &
            fice_file%data_update, fice_file%data_type)

    !---------------------------------------------------------------------------
    !   next, read piston velocity file
    !---------------------------------------------------------------------------

       allocate(xkw_file%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
       in_file = construct_file(xkw_file%input%file_fmt,                 &
                         full_name=trim(xkw_file%input%filename), &
                         record_length = rec_type_dbl,                &
                         recl_words=nx_global*ny_global)
       call data_set(in_file,'open_read')

       i_dim = construct_io_dim('i',nx_global)
       j_dim = construct_io_dim('j',ny_global)
       month_dim = construct_io_dim('month',12)

       io_tracer = &
           construct_io_field(trim(xkw_file%input%file_varname), &
           dim1=i_dim, dim2=j_dim, dim3=month_dim,               &
           field_loc = field_loc_center,                         &
           field_type = field_type_scalar,                       &
           d3d_array=TEMP_DATA)

       call data_set(in_file,'define',io_tracer)
       call data_set(in_file,'read'  ,io_tracer)
       call destroy_io_field(io_tracer)

!untested       !$OMP PARALLEL DO PRIVATE(iblock, n)
       do iblock=1,nblocks_clinic
       do n=1,12
          xkw_file%DATA(:,:,iblock,1,n) = &
             TEMP_DATA(:,:,n,iblock)*xkw_file%input%scale_factor
          where (.not. LAND_MASK(:,:,iblock))
             xkw_file%DATA(:,:,iblock,1,n) = c0
          end where
       end do
       end do
!untested       !$OMP END PARALLEL DO

       call data_set(in_file,'close')
       call destroy_file(in_file)

       if (my_task == master_task) then
          write(stdout,blank_fmt)
          write(stdout,'(a12,a)') ' file read: ', &
             trim(xkw_file%input%filename)
       endif

       call find_forcing_times(xkw_file%data_time, &
            xkw_file%data_inc, xkw_file%interp_type, &
            xkw_file%data_next, xkw_file%data_time_min_loc, &
            xkw_file%data_update, xkw_file%data_type)

    !---------------------------------------------------------------------------
    !   last, read atmospheric pressure file
    !---------------------------------------------------------------------------

       allocate(ap_file%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
       in_file = construct_file(ap_file%input%file_fmt,                &
                         full_name=trim(ap_file%input%filename),       &
                         record_length = rec_type_dbl,                 &
                         recl_words=nx_global*ny_global)
       call data_set(in_file,'open_read')

       i_dim = construct_io_dim('i',nx_global)
       j_dim = construct_io_dim('j',ny_global)
       month_dim = construct_io_dim('month',12)

       io_tracer = &
           construct_io_field(trim(ap_file%input%file_varname), &
           dim1=i_dim, dim2=j_dim, dim3=month_dim,              &
           field_loc = field_loc_center,                        &
           field_type = field_type_scalar,                      &
           d3d_array=TEMP_DATA)

       call data_set(in_file,'define',io_tracer)
       call data_set(in_file,'read'  ,io_tracer)
       call destroy_io_field(io_tracer)

!untested       !$OMP PARALLEL DO PRIVATE(iblock, n)
       do iblock=1,nblocks_clinic
       do n=1,12
          ap_file%DATA(:,:,iblock,1,n) = &
             TEMP_DATA(:,:,n,iblock)*ap_file%input%scale_factor
          where (.not. LAND_MASK(:,:,iblock))
             ap_file%DATA(:,:,iblock,1,n) = c0
          end where
       end do
       end do
!untested       !$OMP END PARALLEL DO

       call data_set(in_file,'close')
       call destroy_file(in_file)

       if (my_task == master_task) then
          write(stdout,blank_fmt)
          write(stdout,'(a12,a)') ' file read: ', &
             trim(ap_file%input%filename)
       endif

       call find_forcing_times(ap_file%data_time, &
            ap_file%data_inc, ap_file%interp_type, &
            ap_file%data_next, ap_file%data_time_min_loc, &
            ap_file%data_update, ap_file%data_type)

       deallocate(TEMP_DATA)

     case ('bulk')

       if (my_task == master_task) then
          write(stdout,*)  &
             ' Using fields from bulk forcing for calculating CFC11 flux'
       endif

     case default
      call exit_POP(sigAbort, &
                    'cfc11_init_sflux: Unknown value for cfc11_formulation')

   end select

    !---------------------------------------------------------------------------
    !   now read latitudinally dependent weights file
    !---------------------------------------------------------------------------

       allocate(CFC11_INTERP_WEIGHT(nx_block,ny_block,max_blocks_clinic))
       in_file = construct_file(cfc11_interp_weight_file_fmt,                 &
                         full_name=trim(cfc11_interp_weight_file), &
                         record_length = rec_type_dbl,                &
                         recl_words=nx_global*ny_global)
       call data_set(in_file,'open_read')

       i_dim = construct_io_dim('i',nx_global)
       j_dim = construct_io_dim('j',ny_global)

       io_tracer = &
           construct_io_field(trim(ap_file%input%file_varname), &
           dim1=i_dim, dim2=j_dim,                              &
           field_loc = field_loc_center,                        &
           field_type = field_type_scalar,                      &
           d2d_array=CFC11_INTERP_WEIGHT)

       call data_set(in_file,'define',io_tracer)
       call data_set(in_file,'read'  ,io_tracer)
       call destroy_io_field(io_tracer)
       call data_set(in_file,'close')
       call destroy_file(in_file)

       if (my_task == master_task) then
          write(stdout,blank_fmt)
          write(stdout,'(a12,a)') ' file read: ', &
             trim(cfc11_interp_weight_file)
       endif

    !---------------------------------------------------------------------------
    !   now read in hemispherically averaged CFC11 time series file
    !   first, read until end of file to determine how big to allocate
    !      the time series array
    !---------------------------------------------------------------------------

      call get_unit(nu)
      if (my_task == master_task) then
        write(stdout,*) 'Trying to read CFC11 time series from file:',  &
                           trim(cfc11_time_series_file)
        open(nu, file=cfc11_time_series_file, status='old', iostat= io_error)
        if (io_error /= 0) then
           io_error = -1
        else
           io_error =  1
        endif
        allocate( cfc11_time_series(1,3) )   ! do this temporarily
        num_cfc11_years = 0
        do while (io_error >= 0)
           read(nu,*,iostat = io_error) ( cfc11_time_series(1,nc), nc = 1, 3 )
           num_cfc11_years = num_cfc11_years + 1
        end do
        num_cfc11_years = num_cfc11_years - 1
        close (nu)
        write(stdout,*)num_cfc11_years,' years to be read from file: ',  &
                       trim(cfc11_time_series_file)
        deallocate( cfc11_time_series )
      endif

    !---------------------------------------------------------------------------
    !   check that num_cfc11_years is > 0
    !---------------------------------------------------------------------------

    call broadcast_scalar(num_cfc11_years, master_task)
    if (num_cfc11_years == 0) then
       call exit_POP(sigAbort,'cfc11_init_sflux: num_cfc11_years = 0')
    end if

    !---------------------------------------------------------------------------
    !   allocate time series array on all processors, then read in on
    !      master_task, then broadcast to all
    !   time series file has 3 columns:  year, Northern Hemisphere, Southern Hem
    !---------------------------------------------------------------------------

      allocate( cfc11_time_series(num_cfc11_years,3) )

      if (my_task == master_task) then
        open(nu, file=cfc11_time_series_file, status='old', iostat= io_error)
        do ny = 1, num_cfc11_years
           read(nu,*,iostat = io_error) ( cfc11_time_series(ny,nc), nc = 1, 3 )
        end do
        close (nu)
        write(stdout,*)' done reading file: ', trim(cfc11_time_series_file)
      endif
      call release_unit(nu)

      call broadcast_array (cfc11_time_series, master_task)

!-----------------------------------------------------------------------
!EOC

 end subroutine cfc11_init_sflux

!***********************************************************************
!BOP
! !IROUTINE: cfc11_set_sflux
! !INTERFACE:

 subroutine cfc11_set_sflux(WIND_VEL,IFRAC,PRESS,SST,SSS,SURF_VALS,  &
                              STF_MODULE)

! !DESCRIPTION:
!  Compute CFC11 surface flux and store related tavg fields for
!  subsequent accumulating.

! !REVISION HISTORY:
!  same as module

! !USES:

   use constants, only: field_loc_center, field_type_scalar, &
       mpercm, eps, p5, cmperm
   use time_management, only: thour00, check_time_flag, iyear,  &
       seconds_this_year, seconds_in_year, tday
   use broadcast, only: broadcast_scalar
   use forcing_tools, only: update_forcing_data, interpolate_forcing

!-----------------------------------------------------------------------
!   NOTE that variables are a mish-mash of model/non-model units
!   WIND_VEL should be m/s
!   XKW should be cm/s
!   salinity should be psu
!   pressure should be dyne/cm**2 (NOT atmospheres or Pascals)
!-----------------------------------------------------------------------

! !INPUT PARAMETERS:

   real(r8), dimension(nx_block,ny_block,2,max_blocks_clinic), intent(in) :: &
         WIND_VEL ! surface wind velocity (m/s)

   real(r8), dimension(nx_block,ny_block,max_blocks_clinic), intent(in) :: &
         IFRAC, & ! sea ice fraction (non-dimensional)
         PRESS, & ! sea level atmospheric pressure (dyne/cm**2)
         SST,   & ! sea surface temperature (C)
         SSS      ! sea surface salinity (psu)

   real(r8), dimension(nx_block,ny_block,cfc11_tracer_cnt,max_blocks_clinic), &
         intent(in) :: SURF_VALS ! module tracers

! !OUTPUT PARAMETERS:

   real(r8), dimension(nx_block,ny_block,cfc11_tracer_cnt,max_blocks_clinic), &
         intent(inout) :: STF_MODULE

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

    character(*), parameter :: subname = 'cfc11_mod:cfc11_set_sflux'

    integer(int_kind) :: &
         j, iblock,    & ! loop indices
         cfc11_first_year, cfc11_last_year  !  first, last year in cfc11 time series

    real(r8) :: &  !  used for linear interpolation of time series values
          pcfc11_north, pcfc11_south, del_year

    real(r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
         FICE_USED,    & ! used ice fraction (non-dimensional)
         XKW,          & ! part of piston velocity (cm/s)
         AP_USED         ! used atm pressure (converted from dyne/cm**2 to atm)

    real(r8), dimension(nx_block,ny_block) :: &
         CFC11_loc,    & ! local copy of surface cfc11, forced to be non-negative
         XKW_ICE,      & ! common portion of piston vel., (1-fice)*xkw (cm/s)
         SCHMIDT_USED, & ! used Schmidt number
         PV,           & ! piston velocity (cm/s)
         CFC11SAT_1atm,& ! calculated CFC11 saturation (mmol/m^3)
         CFC11SAT_USED,& ! used CFC11 saturation (mmol/m^3)
         FLUX            ! tracer flux (nmol/cm^2/s)

   character (char_len) :: &
      tracer_data_label          ! label for what is being updated

   character (char_len), dimension(:), allocatable :: &
      tracer_data_names          ! short names for input data fields

   integer (int_kind), dimension(:), allocatable :: &
      tracer_bndy_loc,          &! location and field type for ghost
      tracer_bndy_type           !    cell updates

   real (r8), dimension(:,:,:,:), allocatable :: &
      TEMP_DATA           ! temp array for reading Tracer data

    real(r8), dimension(nx_block,ny_block) :: &
         WORK1, WORK2 ! temporaries for averages

    !---------------------------------------------------------------------------
    !   local parameters
    !---------------------------------------------------------------------------

    real(r8), parameter :: &
         xkw_coeff = 8.6e-9_r8       ! xkw_coeff = 0.31 cm/hr s^2/m^2 in (s/cm)

    !---------------------------------------------------------------------------
    !---------------------------------------------------------------------------

    STF_MODULE = c0

    !---------------------------------------------------------------------------
    !   Interpolate gas flux forcing data if necessary
    !---------------------------------------------------------------------------

    select case (cfc11_formulation)

    case ('ocmip')

       if (thour00 >= fice_file%data_update) then
          allocate(tracer_data_names(1), &
                   tracer_bndy_loc  (1), &
                   tracer_bndy_type (1))
          tracer_data_names = fice_file%input%file_varname
          tracer_bndy_loc   = field_loc_center
          tracer_bndy_type  = field_type_scalar
          tracer_data_label = 'Ice Fraction'
          call update_forcing_data(          fice_file%data_time,      &
               fice_file%data_time_min_loc,  fice_file%interp_type,    &
               fice_file%data_next,          fice_file%data_update,    &
               fice_file%data_type,          fice_file%data_inc,       &
               fice_file%DATA(:,:,:,:,1:12), fice_file%data_renorm,    &
               tracer_data_label,            tracer_data_names,        &
               tracer_bndy_loc,              tracer_bndy_type,         &
               fice_file%filename,           fice_file%input%file_fmt)
          deallocate(tracer_data_names, &
                     tracer_bndy_loc  , &
                     tracer_bndy_type )
       end if
       allocate(TEMP_DATA(nx_block,ny_block,max_blocks_clinic,1))
       call interpolate_forcing(TEMP_DATA, &
            fice_file%DATA(:,:,:,:,1:12), &
            fice_file%data_time,         fice_file%interp_type, &
            fice_file%data_time_min_loc, fice_file%interp_freq, &
            fice_file%interp_inc,        fice_file%interp_next, &
            fice_file%interp_last,       0)
       FICE_USED = TEMP_DATA(:,:,:,1)
       deallocate(TEMP_DATA)

       if (thour00 >= xkw_file%data_update) then
          allocate(tracer_data_names(1), &
                   tracer_bndy_loc  (1), &
                   tracer_bndy_type (1))
          tracer_data_names = xkw_file%input%file_varname
          tracer_bndy_loc   = field_loc_center
          tracer_bndy_type  = field_type_scalar
          tracer_data_label = 'Piston Velocity'
          call update_forcing_data(         xkw_file%data_time,      &
               xkw_file%data_time_min_loc,  xkw_file%interp_type,    &
               xkw_file%data_next,          xkw_file%data_update,    &
               xkw_file%data_type,          xkw_file%data_inc,       &
               xkw_file%DATA(:,:,:,:,1:12), xkw_file%data_renorm,    &
               tracer_data_label,           tracer_data_names,        &
               tracer_bndy_loc,             tracer_bndy_type,         &
               xkw_file%filename,           xkw_file%input%file_fmt)
          deallocate(tracer_data_names, &
                     tracer_bndy_loc  , &
                     tracer_bndy_type )
       end if
       allocate(TEMP_DATA(nx_block,ny_block,max_blocks_clinic,1))
       call interpolate_forcing(TEMP_DATA,     &
            xkw_file%DATA(:,:,:,:,1:12), &
            xkw_file%data_time,         xkw_file%interp_type, &
            xkw_file%data_time_min_loc, xkw_file%interp_freq, &
            xkw_file%interp_inc,        xkw_file%interp_next, &
            xkw_file%interp_last,       0)
       XKW = TEMP_DATA(:,:,:,1)
       deallocate(TEMP_DATA)

       if (thour00 >= ap_file%data_update) then
          allocate(tracer_data_names(1), &
                   tracer_bndy_loc  (1), &
                   tracer_bndy_type (1))
          tracer_data_names = ap_file%input%file_varname
          tracer_bndy_loc   = field_loc_center
          tracer_bndy_type  = field_type_scalar
          tracer_data_label = 'Atmospheric Pressure'
          call update_forcing_data(        ap_file%data_time,      &
               ap_file%data_time_min_loc,  ap_file%interp_type,    &
               ap_file%data_next,          ap_file%data_update,    &
               ap_file%data_type,          ap_file%data_inc,       &
               ap_file%DATA(:,:,:,:,1:12), ap_file%data_renorm,    &
               tracer_data_label,          tracer_data_names,        &
               tracer_bndy_loc,            tracer_bndy_type,         &
               ap_file%filename,           ap_file%input%file_fmt)
          deallocate(tracer_data_names, &
                     tracer_bndy_loc  , &
                     tracer_bndy_type )
       end if
       allocate(TEMP_DATA(nx_block,ny_block,max_blocks_clinic,1))
       call interpolate_forcing(TEMP_DATA, &
            ap_file%DATA(:,:,:,:,1:12), &
            ap_file%data_time,         ap_file%interp_type, &
            ap_file%data_time_min_loc, ap_file%interp_freq, &
            ap_file%interp_inc,        ap_file%interp_next, &
            ap_file%interp_last,       0)
       AP_USED = TEMP_DATA(:,:,:,1)
       deallocate(TEMP_DATA)

   case ('bulk')

!untested      !$OMP PARALLEL DO PRIVATE(iblock)
      do iblock = 1, nblocks_clinic

         FICE_USED(:,:,iblock) = IFRAC(:,:,iblock)
         XKW(:,:,iblock) = cmperm*cmperm * xkw_coeff*  &  ! wind is m/s, xkw is cm/s
            (WIND_VEL(:,:,1,iblock)**2 + WIND_VEL(:,:,2,iblock)**2)
         AP_USED(:,:,iblock) = PRESS(:,:,iblock)
      enddo
!untested      !$OMP END PARALLEL DO

   end select

       !------------------------------------------------------------------------
       !   assume PRESS is in cgs units (dyne/cm**2) since that is what is
       !     required for pressure forcing in barotropic
       !   want units to be atmospheres
       !   convertion from dyne/cm**2 to Pascals is P(mks) = P(cgs)/10.
       !   convertion from Pascals to atm is P(atm) = P(Pa)/101.325e+3_r8
       !
       !   Set bad AP values to 1. This is necessary for runs restarting off
       !   a run in which the flux coupler didnot restart on AP correctly.
       !------------------------------------------------------------------------

!untested   !$OMP PARALLEL DO PRIVATE(iblock)
   do iblock = 1, nblocks_clinic

      AP_USED(:,:,iblock) = AP_USED(:,:,iblock) / 101.325e+4_r8

      AP_USED(:,:,iblock) = merge( c1, AP_USED(:,:,iblock),   &
                                   (AP_USED(:,:,iblock) > 1.5_r8 .or.   &
                                    AP_USED(:,:,iblock) < 0.5_r8) )
   enddo
!untested   !$OMP END PARALLEL DO

    !---------------------------------------------------------------------------
    !   make sure current model year is within range of cfc11 time series
    !---------------------------------------------------------------------------

    cfc11_first_year = int(cfc11_time_series(1,1) + eps)
    cfc11_last_year  = int(cfc11_time_series(num_cfc11_years,1) + eps)
    if (iyear < cfc11_first_year .or. iyear > cfc11_last_year ) then
       call exit_POP(sigAbort,  &
                'model year out of range of CFC11 years')
    end if

    !---------------------------------------------------------------------------
    !   use linear interpolation of time series to current time
    !---------------------------------------------------------------------------

    del_year = seconds_this_year/seconds_in_year
    j = iyear - cfc11_first_year + 1

    if (iyear == cfc11_first_year .and. del_year <= p5) then
       pcfc11_north = cfc11_time_series(1,2)
       pcfc11_south = cfc11_time_series(1,3)
    elseif (iyear == cfc11_last_year .and. del_year >= p5) then
       pcfc11_north = cfc11_time_series(num_cfc11_years,2)
       pcfc11_south = cfc11_time_series(num_cfc11_years,3)
    else
       if (del_year >= p5) then
          pcfc11_north =  &
             (1.5_r8 - del_year)*cfc11_time_series(j    ,2) +  &
             (del_year - p5    )*cfc11_time_series(j + 1,2)
          pcfc11_south =  &
             (1.5_r8 - del_year)*cfc11_time_series(j    ,3) +  &
             (del_year - p5    )*cfc11_time_series(j + 1,3)
       else
          pcfc11_north =  &
             (p5 + del_year)*cfc11_time_series(j    ,2) +  &
             (p5 - del_year)*cfc11_time_series(j - 1,2)
          pcfc11_south =  &
             (p5 + del_year)*cfc11_time_series(j    ,3) +  &
             (p5 - del_year)*cfc11_time_series(j - 1,3)
       endif
    endif

!maltrud debug
!   if (my_task == master_task) write(stdout,*)  &
!      tday, pcfc11_north, pcfc11_south, 'pCFC11_NS'

    pCFC11 = CFC11_INTERP_WEIGHT*pcfc11_north +   &
             (c1 - CFC11_INTERP_WEIGHT)*pcfc11_south

    !---------------------------------------------------------------------------
    !   Compute XKW_ICE. XKW is zero over land, so XKW_ICE is too.
    !---------------------------------------------------------------------------

!untested    !$OMP PARALLEL DO PRIVATE(iblock,XKW_ICE,SCHMIDT_USED,CFC11SAT_1atm, &
!untested    !$OMP                     CFC11SAT_USED,PV,FLUX)

    do iblock = 1, nblocks_clinic

       XKW_ICE = XKW(:,:,iblock)
       where (FICE_USED(:,:,iblock) > 0.2_r8   &
        .and. FICE_USED(:,:,iblock) < 0.9999_r8)
          XKW_ICE = (c1 - FICE_USED(:,:,iblock)) * XKW_ICE
       end where
       where (FICE_USED(:,:,iblock) >= 0.9999_r8)
          XKW_ICE = c0
       end where

    !---------------------------------------------------------------------------
    !   compute CFC11 flux
    !---------------------------------------------------------------------------

       SCHMIDT_USED = SCHMIDT_CFC(SST(:,:,iblock), LAND_MASK(:,:,iblock), 11)

       CFC11SAT_1atm = pCFC11(:,:,iblock)*   &
          SOLUBILITY_CFC(SST(:,:,iblock),  &
                         SSS(:,:,iblock),   &
                         LAND_MASK(:,:,iblock), 11)

       where (LAND_MASK(:,:,iblock))
          PV = XKW_ICE * SQRT(660.0_r8 / SCHMIDT_USED)
          CFC11SAT_USED = AP_USED(:,:,iblock) * CFC11SAT_1atm
!         CFC11_loc = max(c0, SURF_VALS(:,:,cfc11_ind,iblock))
!         FLUX = PV * (CFC11SAT_USED - CFC11_loc)
          FLUX = PV * (CFC11SAT_USED - SURF_VALS(:,:,cfc11_ind,iblock))
          STF_MODULE(:,:,cfc11_ind,iblock) = FLUX
       end where

       CFC11_SFLUX_TAVG(:,:,1,iblock) = FICE_USED(:,:,iblock)
       CFC11_SFLUX_TAVG(:,:,2,iblock) = XKW_ICE(:,:)
       CFC11_SFLUX_TAVG(:,:,3,iblock) = AP_USED(:,:,iblock)
       CFC11_SFLUX_TAVG(:,:,4,iblock) = pCFC11(:,:,iblock)
       CFC11_SFLUX_TAVG(:,:,5,iblock) = PV * mpercm
       CFC11_SFLUX_TAVG(:,:,6,iblock) = SCHMIDT_USED
       CFC11_SFLUX_TAVG(:,:,7,iblock) = CFC11SAT_USED
!      CFC11_SFLUX_TAVG(:,:,7,iblock) = CFC11_INTERP_WEIGHT(:,:,iblock)

    enddo

!untested    !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!EOC

 end subroutine cfc11_set_sflux

!***********************************************************************
!BOP
! !IROUTINE: cfc11_tavg_forcing
! !INTERFACE:

 subroutine cfc11_tavg_forcing

! !DESCRIPTION:
!  Make accumulation calls for forcing related tavg fields. This is
!  necessary because the forcing routines are called before tavg flags
!  are set.

! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      iblock              ! block loop index

!-----------------------------------------------------------------------

!untested   !$OMP PARALLEL DO PRIVATE(iblock)

   do iblock = 1,nblocks_clinic
      call accumulate_tavg_field(CFC11_SFLUX_TAVG(:,:,1,iblock),tavg_FICE,iblock,1)
      call accumulate_tavg_field(CFC11_SFLUX_TAVG(:,:,2,iblock),tavg_XKW,iblock,1)
      call accumulate_tavg_field(CFC11_SFLUX_TAVG(:,:,3,iblock),tavg_ATM_PRESS,iblock,1)
      call accumulate_tavg_field(CFC11_SFLUX_TAVG(:,:,4,iblock),tavg_pCFC11,iblock,1)
      call accumulate_tavg_field(CFC11_SFLUX_TAVG(:,:,5,iblock),tavg_PV,iblock,1)
      call accumulate_tavg_field(CFC11_SFLUX_TAVG(:,:,6,iblock),tavg_SCHMIDT_CFC11,iblock,1)
      call accumulate_tavg_field(CFC11_SFLUX_TAVG(:,:,7,iblock),tavg_CFC11SAT,iblock,1)
   end do

!untested   !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!EOC

 end subroutine cfc11_tavg_forcing

!***********************************************************************
!BOP
! !IROUTINE: SCHMIDT_CFC
! !INTERFACE:

 function SCHMIDT_CFC(PT_2D,LAND_MASK,kn)

! !DESCRIPTION:
!  CFC 11 and 12 schmidt number
!  as a fonction of temperature.
!
!  ref: Zheng et al (1998), JGR, vol 103,No C1
!
!  Original author: J-C Dutay - LSCE
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   real(r8), dimension(nx_block,ny_block) :: &
      PT_2D              ! temperature (degree Celsius)

   logical (log_kind), dimension(nx_block,ny_block) :: LAND_MASK

   integer(int_kind) ::  &
      kn                 ! 11 = CFC-11, 12 = CFC-12

! !OUTPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block) :: &
      SCHMIDT_CFC        ! returned value, non-dimensional

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   real (r8), dimension(11:12) :: a1, a2, a3, a4

!-----------------------------------------------------------------------
!   coefficients with t in degree Celsius
!-----------------------------------------------------------------------

   a1(11) = 3501.8_r8
   a2(11) = -210.31_r8
   a3(11) =    6.1851_r8
   a4(11) =   -0.07513_r8

   a1(12) = 3845.4_r8
   a2(12) = -228.95_r8
   a3(12) =    6.1908_r8
   a4(12) =   -0.067430_r8

   where (LAND_MASK)
      SCHMIDT_CFC = a1(kn) + a2(kn) * PT_2D + a3(kn) *PT_2D*PT_2D  &
                     + a4(kn) *PT_2D*PT_2D*PT_2D
   elsewhere
      SCHMIDT_CFC = c0
   endwhere

!-----------------------------------------------------------------------
!EOC

 end function SCHMIDT_CFC

!***********************************************************************
!BOP
! !IROUTINE: SOLUBILITY_CFC
! !INTERFACE:

 function SOLUBILITY_CFC(PT,PS,LAND_MASK,kn)

! !DESCRIPTION:
!  CFC 11 and 12 Solubilities in seawater
!  ref: Warner & Weiss (1985) , Deep Sea Research, vol32
!
!  Original author: J-C Dutay - LSCE
!
! !REVISION HISTORY:
!  same as module

! !USES:

   use constants, only: T0_Kelvin, c1000

! !INPUT PARAMETERS:

   real(r8), dimension(nx_block,ny_block) :: &
      PT,              & ! temperature (degree Celsius)
      PS                 ! salinity    (o/oo) (multiply by 1000 to get psu)

   logical (log_kind), dimension(nx_block,ny_block) :: LAND_MASK

   integer(int_kind) ::  &
      kn                 ! 11 = CFC-11, 12 = CFC-12

! !OUTPUT PARAMETERS:

   real(r8), dimension(nx_block,ny_block) :: &
      SOLUBILITY_CFC     ! returned value, in mol/m3/pptv
                         ! 1 pptv = 1 part per trillion = 10^-12 atm = 1 picoatm

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   real(r8), dimension(nx_block,ny_block) :: WORK

   real(r8), dimension(11:12) ::  &
      a1, a2, a3, a4, b1, b2, b3

!-----------------------------------------------------------------------
!  for CFC 11
!-----------------------------------------------------------------------

   a1 ( 11) = -229.9261_r8
   a2 ( 11) =  319.6552_r8
   a3 ( 11) =  119.4471_r8
   a4 ( 11) =   -1.39165_r8
   b1 ( 11) =   -0.142382_r8
   b2 ( 11) =    0.091459_r8
   b3 ( 11) =   -0.0157274_r8

!-----------------------------------------------------------------------
!  for CFC/12
!-----------------------------------------------------------------------

   a1 ( 12) = -218.0971_r8
   a2 ( 12) =  298.9702_r8
   a3 ( 12) =  113.8049_r8
   a4 ( 12) =   -1.39165_r8
   b1 ( 12) =   -0.143566_r8
   b2 ( 12) =    0.091015_r8
   b3 ( 12) =   -0.0153924_r8

   WORK = merge( ((PT + T0_Kelvin)* 0.01_r8), c1, LAND_MASK)

!-----------------------------------------------------------------------
!  coefficient for solubility in  mol/l/atm
!-----------------------------------------------------------------------

   where (LAND_MASK)
      SOLUBILITY_CFC  &
         = exp ( a1 ( kn)                &
               + a2 ( kn)/ WORK          &
               + a3 ( kn)* log ( WORK )  &
               + a4 ( kn)* WORK * WORK   &
               + PS*( ( b3( kn)*WORK + b2( kn) )*WORK + b1(kn) ) )
   elsewhere
      SOLUBILITY_CFC = c0
   endwhere

!-----------------------------------------------------------------------
!  conversion from mol/(l * atm) to mol/(m^3 * atm) to mol/(m3 * pptv)
!-----------------------------------------------------------------------

   SOLUBILITY_CFC = c1000 * SOLUBILITY_CFC
   SOLUBILITY_CFC = 1.0e-12_r8 * SOLUBILITY_CFC

!-----------------------------------------------------------------------
!EOC

 end function SOLUBILITY_CFC

!***********************************************************************

end module cfc11_mod

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
